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1  Major  research  accomplishments 

1.1  A  new  shock-fitting  algorithm  for  reactive  Euler  equations 

We  have  developed  a  new  shock-fitting  algorithm  that  is  based  on  a  coordinate  transformation  to  a 
shock  attached  frame.  We  derived  a  new  evolution  system  for  the  shock  dynamics  from  the  analysis 
of  the  Rankine-Hugoniot  conditions  and  used  it  coupled  with  the  Euler  equations  behind  the  shock 
to  solve  for  the  dynamics  of  detonation  with  high  accuracy.  Our  algorithm  eliminates  numerical 
errors  due  to  shock  capturing  (i.e.  Gibbs  oscillations  and  shock  smearing)  at  the  lead  shock  and 
thus  makes  the  order  of  accuracy  of  the  method  controlled  purely  by  the  discretization  algorithm  in 
the  post-shock  flow.  Our  present  implementation  is  a  high-order  WENO  shock-capturing  method 
for  the  entire  system.  Most  of  the  calculations  reported  in  our  main  paper  [17]  use  a  third  order 
spatial  and  temporal  discretizations  of  the  whole  system  and  show  overall  third  order  accuracy. 

The  main  idea  of  the  proposed  numerical  approach  is  based  on  a  coordinate  transfromation 
of  the  reactive  Euler  equations  to  a  frame  attached  to  the  lead  detonation  shock  with  subsequent 
solution  of  the  postshock  flow  as  well  as  the  evolution  of  the  lead  shock.  A  key  advantage  of  the 
method  is  that  the  coordinate  transformation  is  exactly  the  same  as  in  the  theories  of  detonation 
linear  stability  but  without  linearization.  This  allows  for  a  simple  and  direct  calculation  of  linear 
stability  properties  and  their  comparison  with  theoretical  predictions. 

The  reactive  Euler  equations  of  inviscid  gas  dynamics  ([17,  6]): 


v,  4-  u  •  Vv  —  vV  •  u  =  0, 

(1) 

u,  +u- Vu  +  vVp  =  0, 

(2) 

pt  +u-  Vp  +  pc^V -u  =  pc2O(0, 

(3) 

Af+u-VA  =  (0, 

(4) 

are  coupled  with  the  Rankine-Hugoniot  conditions 

M  —  PaDn  —  PsUnsi 

(5) 

P  —  Pa  +  M2Va  =  ps  +M2Vs, 

(6) 

Dxi  =  Uxis,  Dx 2  —  Uz2si 

(7) 

H  =  ea  +  pava  +  ^ M2vl  =  es+ psvs  + 

(8) 

where  D„=Mva  is  the  speed  of  the  shock  normal  to  itself,  Dx\  and  Dti  are  the  tangential  velocity 
components  of  the  incoming  flow  for  an  observer  moving  with  the  shock,  and  Un,  Ux\,  U?i  are 
the  components  of  the  post-shock  particle  speed  relative  to  the  shock.  Subscripts  a  and  s  identify 
the  ambient  state  immediately  ahead  of  the  shock  and  the  state  immediately  behind  the  shock, 
respectively.  The  shock  state  is  completely  determined  by  the  normal  mass  flux  M  and  the  shock 
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slope  (which  enters  the  definitions  of  the  velocity  components.)  Therefore,  the  shock  dynamics 
will  be  completely  determined  if  we  can  solve  for  the  postshock  flow  as  well  as  the  shock  state, 
given  by  M  and  the  shock  slope. 

If  <i>(x,y,z,t )  =  0  is  the  equation  of  the  shock  surface,  then  the  evolution  of  the  surface  in 
space-time  is  governed  by  the  level-set  equation 


+  D  •  V  O  —  0, 


(9) 


where  D  is  the  surface  velocity  vector.  The  unit  outward  normal  to  the  surface  is  given  by  n  = 
—VO/  |VO|,  and  with  D  =  Dn n  the  normal  mass  flux  across  the  surface  is 


M  = 


(10) 


The  second  equation  needed  for  the  shock  evolution  calculation  comes  from  differentiating  the 
Rankine-Hugoniot  conditions.  We  differentiate  the  post-shock  state  variables,  ps  and  u5,  in  time 
and  relate  their  derivatives  to  the  flow  gradients  at  the  shock  by  means  of  the  Euler  equations.  We 
arrive  at 

d-Ei-psVs-^=  ,[(y- 1  )Qp<o- pc2V •  u +pU •  (U ■  V) u]  |„  (11) 

where  u  —  D  =  U  is  the  particle  velocity  relative  to  the  shock. 

Equation  (11)  together  with  (10)  can  be  solved  for  M  and  O  provided  the  gradient  of  the  post¬ 
shock  velocity  vector  is  known.  Once  M  and  O  are  known,  all  the  shock-state  variables  can  be 
calculated  through  the  Rankine-Hugoniot  conditions.  In  principle,  this  idea  can  be  applied  in  any 
geometry  thus  forming  a  foundation  for  a  generalized  DSD  method,  but  so  far  the  method  has  been 
implemented  only  in  a  channel  geometry  and  a  one-dimensional  spherical/cylindrical  geometries. 

In  channel  geometry,  the  Euler  equations  with  a  one-step  chemical  reaction  are 


dY  dF(Y)  <9G(Y) 
9t  +  dx  +  dy 


S(Y), 


(12) 


where  the  vectors  representing  the  state  variables  Y,  x-fluxes  F,  y-fluxes  G,  and  the  source  term  S 
are 


p 
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Figure  1:  Schematics  of  the  mapping  and  coordinate  transformation  to  the  shock  attached  frame: 
(a)  -  the  actual  shock  shape  in  the  physical  plane;  (b)  -  the  shock  in  the  mapped  computational 
domain. 


In  an  ideal  gas,  the  total  energy  et  is  the  sum  of  the  internal,  chemical,  and  kinetic  energies 
as  e,  =  —  XQ+  .  We  also  take  one-step  irreversible  Arrhenius  reaction  rate  law,  co  = 

k(l  —  A)exp  (-^). 

After  collecting  all  the  terms  in  (1 1)  and  (10),  we  derive  the  following  shock  evolution  system 


where 

'  M 


and 

s  =  s(M,  y/,q)  =  —  (Rs  ~AS),  (16) 

Ao 

A0  =  -i_Mva(3  +  ^,  Rs  =  (Y-1)QPs<Os,  (17) 

A,  =  ps  [(c?  -  U2)  qn  +  (c2  -  V2)  {q22  +  yq2l)  -  USVS  (qn+qil  +  Wll)]  •  (18) 

Here  qy  =  dui/Bxj  are  the  components  of  the  velocity-gradient  tensor  evaluated  at  the  shock 
((jcj  ,x2)  =  (x,y)).  Importantly,  (14)  is  exact.  In  general,  q  is  not  known  in  terms  of  y,  but  we 
approximate  it  based  on  a  previous  time  step  in  a  numerical  algorithm.  Away  from  the  lead  shock. 


y,+fy  =  s, 


(14) 


,  f  = 


-Vay/T+V^M 


,  s  = 


s 

0  ’ 


(15) 
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Density 


Reaction  rate:  t=150 


Figure  2:  Typical  profiles  of  computed  density  and  reaction  rate  in  a  weakly  unstable  detonation 


we  integrate  the  Euler  equations  in  their  conservative  form,  which  is  preserved  under  the  coordinate 
transformation  to  the  shock-attached  frame: 


Yt  + 


F+Mva^l  +  yr2Y  +  y/G 


+  Gy  =  S. 


(19) 


The  algorithm  has  been  implemented  in  parallel  Fortran  90/95  with  MPI.  Most  of  the  calcula¬ 
tions  have  been  carried  out  on  an  88  core  Linux  cluster  at  UIUC  and  took  10-20  hours  to  complete. 
Despite  the  fact  that  the  method  computes  only  the  regions  of  interest  (i.e.  does  not  calculate  the 
unreacted  state  ahead  of  the  shock  or  the  completely  burnt  state  far  behind  it,)  the  computational 
resources  requred  for  well-resolved  simulations  are  still  extensive.  This  is  a  well-known  challenge 
in  detonation  simulations  and  is  a  reflection  of  a  very  complex  dynamics  of  flow  within  the  reaction 
zone.  A  typical  plot  of  the  computed  detonation  structure  is  shown  in  Figure  2. 
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1.2  Validation  of  theoretical  predictions  of  two-dimensional  detonation  in¬ 
stability 

So  far,  we  have  performed  the  first  detailed  testing  of  detonation  stability  theory  in  two  dimensions. 
With  moderate  resolution  of  about  40  points  per  half-reaction  zone,  we  have  been  able  to  compute 
the  growth  rates  and  frequencies  of  instability  in  excellent  agreement  with  theoretical  predictions 
based  on  the  normal-mode  analysis  of  Short  and  Stewart  [14];  our  main  results  appear  in  [17]. 

The  governing  equations  for  the  flow  variables  in  the  shock  attached  frame  and  shock  evolution 
equations  were  discretized  on  a  two-dimensional  finite  volume  grid  and  a  one-dimensional  finite 
volume  grid,  respectively.  Spatial  fluxes  of  the  flow  variables  and  the  flux  in  the  shock  slope 
equation  were  computed  using  a  fifth-order  WENO  scheme  [15].  A  third-order  TVD  Runge-Kutta 
scheme  was  used  for  temporal  discretization.  The  numerical  fluxes  at  the  x-face  nearest  the  shock 
require  special  treatment,  as  we  do  not  want  to  compute  the  flux  interpolants  using  data  from  ahead 
of  the  shock.  In  this  case,  the  interpolant  weights  corresponding  to  the  stencil  that  extends  ahead 
of  the  shock  are  set  to  a  very  small  value  so  that  the  reconstructed  fluxes  are  based  only  on  the 
solution  within  the  computational  domain. 

The  spatial  discretization  scheme  is  shown  schematically  in  Figure  3.  The  figure  illustrates 
that  the  shock  mapping  variables,  M  and  yr,  are  stored  in  a  one-dimensional  grid  along  the  y 
direction.  The  velocity  gradient  tensor  q  in  the  mass  flux  equation  is  computed  using  one-sided 
and  centered  finite-differences  in  the  x  and  y  directions,  respectively.  The  cell-centered  grid  points 
used  to  compute  the  velocity  gradients  in  the  shock-normal  direction  are  indicated  by  filled  black 
circles,  while  the  points  used  to  compute  the  velocity  gradients  in. the  shock-tangential  direction 
are  indicated  by  open  circles. 

A  weakly  unstable  two-dimensional  detonation  with  particularly  simple  unstable  behavior  was 
one  of  the  cases  examined  by  Short  and  Stewart  [14],  and  is  used  as  a  point  of  comparison  here  to 
verify  the  simulation  results.  The  model  parameters  are  y  =  1 .2,  Q  =  0.4,  /  =1.2,  and  E  =  50.  This 
case  is  one-dimensionally  stable  and  two-dimensionally  unstable  over  a  small  range  of  transverse 
disturbance  wavenumbers,  see  Figure  4. 

Simulations  were  performed  in  channels  of  width  2.1  to  7.8  with  solid  wall  boundary  condi¬ 
tions  at  a  resolution  of  30  grid  cells  per  unit  length.  In  each  simulation,  the  longest  admissible 
wavelength  was  perturbed(w  =  2 W).  The  resulting  transverse  disturbance  wavenumbers  cover  the 
entire  unstable  range  for  this  case.  Dispersion  relations  can  be  constructed  from  the  computed 
growth  rates  and  frequencies,  as  shown  in  Figure  4.  As  can  be  seen  from  the  figure,  the  agreement 
with  linear  stability  theory  is  excellent.  Another  example,  now  with  a  more  complex  dispersion 
relation  is  shown  in  Figure  5. 
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Mapping  Variables  (ID) 

Mass  flux,  Shock  slope 


Flow  Variables  (2D) 

Mass,  Momentum,  Energy,  Lambda 


Figure  3:  Spatial  discretization.  Shock  surface  is  gray.  Finite  difference  stencils  for  x  and  y 
velocity  derivatives  at  shock  are  denoted  by  black  and  white  circles,  respectively. 


Figure  4:  Dispersion  relations.  Linear  stability  (solid  line),  simulation  (circles).  Model  parameters: 
y  =  1.2,  2  =  0.4,  /  =  1.2,  E  =  50. 
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Figure  5:  Comparison  of  numerically  computed  spectra  with  linear  stability  results.  Model  param¬ 
eters:  y=  1.2,  g  =  £  =  50,/=  1.8. 


1.3  Extension  of  the  algorithm  to  complex  chemistry 

The  algorithm  can  be  extended  to  include  multiple-step  chemical  mechanisms  by  replacing  in  the 
shock  dynamic  condition  the  energy  release  rate  (product  of  thermicity  and  reaction  rate)  for  the 
single  reaction  with  the  net  energy  release  rate  for  several  reactions  (a  sum  of  all  thermicity  and 
reaction  rate  products.)  This  multi-step  reaction  scheme  was  implemented  in  a  modified  version 
of  SAF2D  code  and  used  to  compute  2D  linear  stability  results  of  a  detonation  model  consisting 
of  an  ideal  EOS  and  a  three-step  chain-braching  reaction  mechanism  [13], 

F-+X,  (Oi  =exp^-£/  > 

F+X-+2X,  ^  =  exp(-£BQ--i^, 

X-+P,  <03  =  1. 

The  shock  dynamic  condition  for  complex  chemistry  is  generalized  as  follows 

^^‘=[pc2Eo-«-pc2V.u+pU(UV)a]s, 

Pi  =  -e)J  (pc2‘p)  , 

while  the  kinematic  condition  is  unchanged.  Required  modifications  for  the  code  are  minimal. 
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Figure  6:  Numerically  computed  spectra  of  linear  instability  for  detonation  with  a  model  three-step 
chain  branching  chemistry  with  £/  =  20,  Eb  =  10,  y  =  1.2,  /  =  1,  Q  =  3,  7)/ 7#  =  3,  Tb/Ts  —  0.9 
[13]. 


A  particular  case  is  shown  in  Figure  6.  These  results  have  been  presented  at  2009  ICDERS 
meeting  [11]  and  are  being  prepared  for  a  journal  publication.  They  are  also  a  part  of  the  doctoral 
thesis  of  B.  Taylor  [16]. 

1.4  Application  to  radially  expanding  detonations 

For  radially  expanding  detonations,  we  need  to  add  a  divergence  term  to  the  thermicity  in  the 
equation  for  the  mass  flux  (16),  so  that  s  becomes  .now  s  =  (Rs  —  psc*usK  —  As)  /Aq,  where  the 
shock  curvature  is  K  =  jjr  with  j  =  1,2  for  cylindrical  and  spherical  detonation,  respectively. 
This  modification  was  implemented  in  SAF1D  and  used  to  verify  DSD  results  by  taking  the  DSD 
solution  as  an  initial  condition  for  the  direct  simulation.  For  example,  the  simulation  results  have 
shown  that  D  ~  0  on  the  D-kappa  curve.  The  scheme  was  also  used  to  study  the  instability  of  DSD 
solutions  that  have  turning  points.  It  was  demostrated  that  the  middle  branch  of  a  D  —  k  curve  with 
turning  points  is  highly  unstable  and,  in  the  case  of  an  ideal  EOS  and  one-step  reaction  mechanism, 
results  in  a  rapid  transition  ("ignition")  to  a  strong  detonation  on  the  upper  branch,  see  Figure  7. 
This  is  the  first  calculation  of  such  a  solution  using  the  DSD  solution  from  the  middle  branch 
as  an  initial  condition.  Ongoing  efforts,  as  part  of  B.  Taylor’s  PhD  thesis,  include  comparisons 
between  the  simulation  results  from  unstable  DSD  solutions  with  linear  stability  theory  of  a  curved 
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p0«-| ,  p0-1 ,  Y-1 A  0=50,  E*50.  v=1  p0-1 ,  Pt«1,Y=1.2. 0=50,  E=20,  v=1 


Figure  7:  Calculations  of  the  dynamics  of  a  curved  detonation  with  SAF1D  and  comparisons  with 
DSD  predictions.  Simulations  were  started  at  the  D—k  curve  and  continued  using  SAF1D.  Both 
ignition  from  the  middle  branch  and  stable  settling  to  the  CJ  detonation  from  the  upper  branch  of 
the.D  —  K  curve  is  observed. 

detonation. 

Our  algorithm  is  ideally  suited  for  direct  calculations  of  detonation  shock  dynamics  using  the 
same  variables  as  those  used  in  DSD  theories  [2,  1,  8,  10],  namely  shock  acceleration  D,  speed  D, 
and  curvature  K  of  the  front.  This  provides  a  useful  tool  for  verification  of  results  of  DSD  theories, 
which  are  otherwise  hard  to  test. 

1.5  Formulation  for  non-ideal  equations  of  state 

Our  shock-attached  frame  formulation  admits  a  generalization  to  a  non-ideal  equation  of  state  as 
well  and  thus  to  realistic  models  of  high  explosives.  Complications  arise  in  the  shock  dynamic 
condition  primarily  because  of  the  need  to  calculate  the  left-hand  side,  which  involves  time  deriva¬ 
tives  of  the  shock  state  variable.  Now  the  left-hand  side  will  be  proportional  to  Mt  with  a  factor 
that  involves  the  terms  dps/dM  and  dus/dM,  which  cannot  be  computed  analytically  for  a  gen¬ 
eral  EOS.  Instead,  the  scheme  implemented  in  SAF1D  numerically  approximates  these  derivatives 
using  a  central  finite  difference  of  solutions  to  the  Rankine-Hugoniot  relations.  This  scheme  was 
tested  against  the  analytical  solution  for  an  ideal  EOS  and  yields  results  that  are  in  very  good  agree¬ 
ment.  Thus  the  only  complications  due  to  non-ideal  EOS  are  of  algebraic  nature  and  require  no 
new  derivations  of  the  shock  evolution  equations.  This  work  is  ongoing  and  is  also  a  part  of  Brian 
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Figure  8:  Computed  modal  content  of  the  detonation  shock  perturbation  showing  the  most  unstable 
dominant  mode  together  with  the  secondary  modes.  The  linear  stability  results,  namely  the  growth 
rate  and  frequency,  are  extracted  from  these  data.  Initially  linear  growth  of  perturbations  evolves 
into  a  non-linear  regime  where  mode  interactions  lead  to  changes  in  the  modal  growth  rates  [17]. 

Taylor’s  PhD  dissertation  [16].  The  non- ideal  EOS  that  is  being  implemented  is  the  wide-ranging 
EOS  of  Wescott,  Stewart,  and  Davis  [18]. 

1.6  A  method  for  direct  calculation  of  detonation  instability  spectra 

One  of  the  most  important  achievements  of  our  research  is  a  new,  direct  method  for  calculation  of 
detonation  linear  stability  properties  for  arbitrary  constitutive  laws  [11].  The  kinematic  condition 
(10)  is  independent  of  the  details  of  chemistry  and  equation  of  state.  At  the  same  time,  as  described 
above,  the  dynamic  condition  can  be  easily  generalized  to  arbitrary  constitutive  laws  [16,  12,  3]. 
Thus  the  method  can  be  used  to  compute  linear  growth  of  instability  of  any  reactive  system  with 
any  chemical  reaction  mechanism.  Although  computational  demands  will  still  be  high,  at  present 
this  is  the  only  available  general  method  for  calculating  detonation  stability  spectra.  It  must  be 
pointed  out  that  the  detonation  stability  theory,  first  formulated  by  Erpenbeck  in  1960s  [4,  5],  is 
still  being  developed  for  general  constitutive  laws  50  years  later.  There  are  major  theoretical  and 
numerical  difficulties  associated  with  the  normal-mode  formulation  of  detonation  stability  theory 
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stemming  from  the  complexity  of  the  far-field  radiation  conditions  and  problems  with  stiff  numer¬ 
ical  ODE  and  root  solvers.  These  problems  are  completely  absent  in  the  direct  numerical  solution 
of  the  stability  problem,  because  our  method  computes  the  evolution  of  small  perturbations  of  the 
full  exact  reactive  Euler  system  employing  the  exact  Rankine-Hugoniot  conditions  and  requiring 
no  radiation  conditions.  In  contrast,  the  normal  mode  analysis  computes  the  solutions  of  the  lin¬ 
earized  Euler  equations  and  requires  knowledge  of  the  complicated  behavior  of  the  solution  in  the 
neighborhood  of  the  frozen  sonic  state.  Our  ongoing  research  efforts  are  directed  to  creating  a 
unified  framework  for  calculation  of  the  reactive  Euler  equations  in  the  shock  attached  frame  in 
1,2,  or  3  dimensions  with  subsequent  postprocessing  that  analyzes  the  simulation  results  to  extract 
the  growth  rates,  frequencies,  and  eigenfunctions  from  the  growth  of  small  perturbations. 

1.7  Development  of  detonation  analogs 

We  have  discovered  that  there  exists  a  number  of  phenomena  that  are  analogous  to  the  Chapman- 
Jouguet  detonation  wave  in  terms  their  structure  consisting  of  a  shock  followed  by  transonic  flow. 
Two  of  such  phenomena,  a  circular  hydraulic  jump  [9]  and  a  traffic  jam  [7],  have  been  investigated 
in  detail.  It  was  shown  that  both  in  traffic  jam  and  the  hydraulic  jump,  the  theory  of  . the  steady 
solution  is  similar  to  the  Zel’dovich-von  Neumann-Doering  theory  of  detonation  [6].  In  both  cases, 
the  steady  structure  is  found  by  solving  the  underlying  hyperbolic  system,  imposing  the  shock 
conditions  and  sonic  (CJ)  conditions.  Unlike  the  detonation  wave,  these  phenomena  are  described 
by  a  system  of  two  partial  differential  equations  and  is  thus  simpler  to  analyze.  Nevertheless, 
going  beyond  the  steady  state  solutions  leads  to  serious  challenges  having  to  do  with  the  behavior 
of  solutions  near  the  sonic  locus.  Our  continuing  efforts  are  directed  towards  understanding  the 
formulation  of  the  initial  boundary  value  problems  for  such  systems  and  the  regularity  of  their 
solutions.  This  appears  to  be  a  major  mathematical  challenge  with  wide-ranging  implications  for 
a  large  number  of  applications  where  self-sustained  shocks  are  possible,  not  only  in  shallow-water 
systems  and  traffic  flows,  but  also  in  certain  astrophysical  flows  (accretion  shocks,)  flows  in  flexible 
tubes,  two-phase  flows,  and  others. 

The  analogy  is  seen  in  the  main  properties  of  the  shallow- water  and  traffic-flow  systems  [9, 7]: 

uf  +F*  =  s. 

In  both  cases,  a  traveling  wave  solution  exists  of  the  form,  u  =  u  (x  —  st),  and  the  resultant  system 
of  ordinary  differential  equations  reduces  to  a  single  differential  equation, 

du  _  fM) 

c2  —  u 2’ 

where  c  is  the  local  speed  of  small  perturbations  and  u  is  the  flow  speed.  This  equation  has  a 
regular  solution  crossing  the  sonic  point  c  =  u  provided  /  —  0  at  the  same  point,  which  is  precisely 
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Figure  9:  Left  -  comparison  between  our  theoretical  prediction  of  the  hydraulic  jump  structure  [9] 
and  experiment.  Right  -  jamiton  profile  in  the  theory  of  traffic  shocks  [7]. 

the  Chapman-Jouguet  condition  of  detonation  theory. 

The  discovery  of  the  analogy  allowed  us  to  develop  the  first  rational  theory  of  a  classical  phe¬ 
nomenon  in  fluid  mechanics,  the  hydraulic  jump  [9],  and  to  find  a  simple  analytical  solution  to 
another  textbook  problem  in  gasdynamics/partial  differential  equations,  namely  the  traffic  jam  [7] 
described  within  the  framework  of  a  continuum  model. 
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Prof.  Ruben  R.  Rosales  (MIT,)  B.  Seibold  (MIT,)  J.-C.  Nave  (MIT,)  and  Prof.  M.  Flynn 
(University  of  Alberta.) 

3  PhD  Thesis 

1.  B.  D.  Taylor,  Stability  of  Steady  and  Quasi-steady  Detonations,  University  of  Illinois  at 
Urbana-Champaign.  Co-advised  with  Prof.  D.S.  Stewart. 
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4  Publications 


1 .  B.  D.  Taylor,  A.  R.  Kasimov,  and  D.  S.  Stewart,  Mode  selection  in  unstable  two-dimensional 
detonations.  Combustion  Theory  and  Modelling,  Vol.  13,  No.  6,  2009,  973-992. 

2.  M.R.  Flynn,  A.  R.  Kasimov,  J.-C.  Nave,  R.R.  Rosales,  B.  Seibold,  Self-sustained  nonlinear 
waves  in  traffic  flow,  Phys.  Rev.  E  79,  0561 13  (2009). 

3.  A.R.  Kasimov,  A  stationary  circular  hydraulic  jump  and  its  gasdynamic  analogue,  Journal 
of  Fluid  Mechanics,  601:189-198,  2008 

4.  A.R.  Kasimov,  Hyperbolic  dynamics  of  weakly  nonlinear  detonation  waves  (in  preparation). 

5.  A.R.  Kasimov,  B.  D.  Taylor,  D.S.  Stewart,  Accurate  direct  numerical  computation  of  deto¬ 
nation  instability,  1CDERS,  Minsk,  2009 

6.  A.R.  Kasimov,  Detonation  analogs,  ICDERS,  Minsk,  2009 

5  Invited  presentations  and  talks 

1 .  A.R.  Kasimov,  Accurate  Simulation  of  Multi-Dimensional  Detonation  Waves  in  a  Shock- 
Attached  Frame,  APS  Division  of  Fluid  Dynamics,  Nov.  2008. 

2.  A.R.  Kasimov,  Detonation:  Its  Theory,  Numerical  Simulation,  and  a  Hydraulic  Analogue, 
(Seminar  at  GALCIT,  Caltech,  Jan  18,  2008). 

3.  A.R.  Kasimov,  Shocks  and  sonic  points,  MIT  Physical  Mathematics  Seminar,  May  13,  2008. 

4.  A.R.  Kasimov,  The  power  of  analogy:  from  traffic  jams  to  black  holes,  MIT  Applied  Math¬ 
ematics  Colloquium,  2009. 

5.  A.R.  Kasimov,  Theory  and  simulation  of  shock  dynamics:  detonations  and  related  phenom¬ 
ena,  MIT  Mechanical  Engineering  Seminar,  2009 

6.  A.R.  Kasimov,  Theory  and  simulation  of  shock  dynamics:  detonations  and  related  phenom¬ 
ena,  NEU  Mechanical  Engineering  Seminar,  2009 

7.  A.R.  Kasimov,  Analysis  and  simulation  of  shock  dynamics:  detonations  and  related  phe¬ 
nomena,  U.  of  Alberta,  2009 

8.  A.R.  Kasimov,  Analysis  and  simulation  of  shock  dynamics:  detonations  and  related  phe¬ 
nomena,  UMass  Dartmouth,  2009 
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6  Conference  talks 


1.  R.  Kasimov,  Kinematics  and  dynamics  of  multi-dimensional  shock  waves  with  application  to 
detonation  simulation,  12th  International  Conference  on  Numerical  Combustion,  CA,  April 
2008 

2.  A.R.  Kasimov,  B.  D.  Taylor,  and  D.  S.  Stewart,  Mode  selection  in  unstable  two-dimensional 
detonations,  12th  International  Conference  on  Numerical  Combustion,  CA,  April  2008 

3.  A.R.  Kasimov,  B.  D.  Taylor,  D.S.  Stewart,  Accurate  direct  numerical  computation  of  deto¬ 
nation  instability,  ICDERS,  Minsk,  2009 

4.  A.R.  Kasimov,  Detonation  analogs,  ICDERS,  Minsk,  2009 

5.  Taylor,  B.  D.,  Kasimov,  A.  R.,  and  Stewart,  D.  S.,  November  2009.  "Axisymmetric  Shock- 
Attached  Frame  Detonation  Simulations."  62nd  Annual  Meeting  of  the  American  Physical 
Society  Division  of  Fluid  Mechanics,  Minneapolis,  Minnesota,  2009 


7  New  discoveries,  inventions,  or  patent  disclosures 

•  A  new  shock  fitting  numerical  algorithm  for  the  simulation  of  multi-dimensional  detonation 
waves 

•  A  new  algorithm  for  calculation  of  detonation  instability  spectra 

•  Analogy  between  detonations,  hydraulic  jumps,  and  traffic  jams 
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